function res=compute_exp(lambda_in,x_in,n,temp_coef,xi,chi,mulambda,alpha,v,mux,rho,phi,ngridpoints)

% user inputs 
% lambda_in: initial value of lambda 
% x_in: initial value of x_t
% n: maturity of the security
% temp_coef: sequence of constants of n elements
% xi: value of a jump
% chi: coefficient on jump indicator in lambda process
% mulambda: drift parameter in lambda process
% alpha: coefficient on lagged lambda value in lambda process
% v: coefficient on lagged x_t value in lambda process
% mux: drfit parameter in x_t process
% rho: coefficient on lagged x_t value in x_t process
% phi: coefficient on jump indicator in x_t process

x_grid = linspace(-1, 1, ngridpoints); % grid of x values
l_grid = linspace(0, 1, ngridpoints); % grid of lambda values

% start of code computing expectation estimate grids
[L,X] = meshgrid(l_grid, x_grid);
current_grid = [];
grids = {};

for i = n:-1:1
    if i == n
        % Calculates range for domain lambda_grid
        current_grid = (L * exp(temp_coef(i)*xi)) + (1 - L);
        grids{n-i+1} = current_grid;
    
    else
        
        % Alternative
        nextL1=mulambda + alpha*L + v*X + chi*xi;
        nextX1=mux + rho*X + phi*xi;
        nextL2=mulambda + alpha*L + v*X;
        nextX2=mux + rho*X;
        
        nextL1(nextL1>1)=1;
        nextL2(nextL2>1)=1;
        nextL1(nextL1<0)=0;
        nextL2(nextL2<0)=0;
        
        f1=interp2(L,X,current_grid,nextL1,nextX1,'linear',NaN);
        f2=interp2(L,X,current_grid,nextL2,nextX2,'linear',NaN);
        
        current_grid = L.*(exp(temp_coef(i)*xi)*f1) +...
                        (1 - L).*f2;
        grids{n-i+1} = current_grid;
        
    end
end

% start of extracting point estimates from grids using initial values of
% lambda_t and x_t
[valL, idL] = min(abs(l_grid - lambda_in));
[valX, idX] = min(abs(x_grid - x_in));
res = nan(n,1);
for i = 1:n
    res(i) = grids{i}(idX,idL);
end